pacman::p_load(sf, tidyverse, tmap, raster, spatstat, maptools, spNetwork, classInt, viridis, arrow, lubridate, dplyr, sfdep)Take-home Exercsie 3 - Prototyping
Installing and Loading R Packages
Importing Geospatial Data
Boundary Layer (Province-Level):
province_sf <- st_read(dsn = "../../data/cambodia/boundary/level1",
layer = "KHM_adm1")Reading layer `KHM_adm1' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\boundary\level1' using driver `ESRI Shapefile'
Simple feature collection with 25 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 102.3355 ymin: 9.91361 xmax: 107.63 ymax: 14.68817
Geodetic CRS: WGS 84st_geometry(province_sf)Geometry set for 25 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 102.3355 ymin: 9.91361 xmax: 107.63 ymax: 14.68817
Geodetic CRS: WGS 84
First 5 geometries:tmap_mode("plot")
tm_shape(province_sf) +
tm_polygons()
Road Layer:
cam_road_sf <- st_read("../../data/cambodia/roads/osm_road_2022_1641440547.gpkg")Reading layer `osm_road_2022_1641440547' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\roads\osm_road_2022_1641440547.gpkg'
using driver `GPKG'
Simple feature collection with 169682 features and 15 fields
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 102.3415 ymin: 10.42264 xmax: 107.6128 ymax: 14.69023
Geodetic CRS: WGS 84unique(cam_road_sf$fclass) [1] "primary" "secondary" "trunk_link" "trunk"
[5] "secondary_link" "primary_link" "unclassified" "residential"
[9] "tertiary" "service" "path" "track"
[13] "steps" "footway" "pedestrian" "track_grade3"
[17] "living_street" "track_grade4" "tertiary_link" "track_grade1"
[21] "cycleway" "track_grade5" "motorway" "motorway_link"
[25] "track_grade2" "bridleway" st_geometry(cam_road_sf)Geometry set for 169682 features
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 102.3415 ymin: 10.42264 xmax: 107.6128 ymax: 14.69023
Geodetic CRS: WGS 84
First 5 geometries:tmap_mode("plot")
tm_shape(cam_road_sf) +
tm_lines()
Healthcare Facilities
Health Centers:
points_healthcenter <- st_read(dsn = "../../data/cambodia/healthcenter",
layer = "healthcenter")Reading layer `healthcenter' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\healthcenter' using driver `ESRI Shapefile'
Simple feature collection with 956 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 213103.9 ymin: 1155396 xmax: 764437 ymax: 1594142
Projected CRS: WGS 84 / UTM zone 48Nhead(points_healthcenter)Simple feature collection with 6 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 276571.9 ymin: 1422232 xmax: 343642 ymax: 1553869
Projected CRS: WGS 84 / UTM zone 48N
PCODE PNAME DCODE DNAME CCODE CNAME VCODE
1 2 Sihanoukville 204 Bavel 20404 Khnach Romeas 2040206
2 2 Battambang 203 Bat Dambang 20301 Tuol Ta aek 2030102
3 2 Battambang 213 Koa Krolor 21301 Thibadey 2130104
4 22 Oddar Mean chey 2203 Chong Kal 220302 Chong kal 22030201
5 1 Banteay Meanchey 102 Mongkol Borei 10201 Banteay Neang 1020103
6 1 Banteay Meanchey 102 Mongkol Borei 10203 Chamnaom 1020305
VNAME ODCODE ODNAME FACILITCOD FACILITNAM
1 Khnach Romeas 0201 Thmar Koul 020117 Khnach Romeas
2 O Takoam II 0204 Bat Dambang 020406 Tuol Ta aek
3 Chhae Balang 0202 Mong Russei 020207 Thibadey
4 Chong Kal 2201 Samraong 220113 Pong ro
5 Banteay Neang 0101 Mongkol Borei 010115 Banteay Neang
6 Chamnaom Kaeut 0101 Mongkol Borei 010112 Chamnaom
COVERNAME geometry
1 <NA> POINT (277999.9 1467055)
2 <NA> POINT (303951.9 1447582)
3 <NA> POINT (312278.9 1422232)
4 <NA> POINT (343642 1553869)
5 Banteay Neang, Commune, 19 Vill POINT (285295.9 1494744)
6 Chamnaom, Commune, 17 Vill POINT (276571.9 1487280)Health Posts:
points_healthpost <- st_read(dsn = "../../data/cambodia/healthpost",
layer = "healthpost")Reading layer `healthpost' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\healthpost' using driver `ESRI Shapefile'
Simple feature collection with 89 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 247772 ymin: 1190286 xmax: 762122.8 ymax: 1591629
Projected CRS: WGS 84 / UTM zone 48Nhead(points_healthpost)Simple feature collection with 6 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 589482.9 ymin: 1326873 xmax: 705301 ymax: 1368201
Projected CRS: WGS 84 / UTM zone 48N
PCODE PNAME DCODE DNAME CCODE CNAME VCODE VNAME
1 10 Kratie 1003 Preaek Prasab 100301 Chambak 10030106 Stueng Thum
2 10 Kratie 1005 Snuol 100501 Khsuem 10050107 Srae Roneam
3 10 Kratie 1005 Snuol 100504 Srae Char 10050404 Mean Chey
4 10 Kratie 1005 Snuol 100505 Svay Chreah 10050508 Rumpuk
5 11 Mondul Kiri 1101 Kaev Seima 110103 Srae Chhuk 11010302 Chak Cha
6 11 Mondul Kiri 1101 Kaev Seima 110104 Srae Khtum 11010401 Ou Am
ODCODE ODNAME FACILITCOD FACILITNAM geometry
1 1001 Chhlong 10010201 Stueng Thum POINT (589482.9 1366181)
2 1002 Kratie 10021501 Srae Roneam POINT (654848 1357465)
3 1002 Kratie 10021401 Doun Meas POINT (644705 1326873)
4 1002 Kratie 10021301 Rumpuk POINT (640357 1352191)
5 1101 Sen Monorom 110108 Srae Chhouk POINT (687516 1368201)
6 1101 Sen Monorom 110108 Ou Am POINT (705301 1339600)Referral Hospitals:
points_referralhospital <- st_read(dsn = "../../data/cambodia/referralhospital",
layer = "hltfacp_referral")Reading layer `hltfacp_referral' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\referralhospital' using driver `ESRI Shapefile'
Simple feature collection with 76 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 216012.9 ymin: 1158852 xmax: 737796 ymax: 1573927
Projected CRS: WGS 84 / UTM zone 48Nhead(points_referralhospital)Simple feature collection with 6 features and 12 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 286126.9 ymin: 1308056 xmax: 628526 ymax: 1542390
Projected CRS: WGS 84 / UTM zone 48N
PCODE PNAME DCODE DNAME CCODE CNAME VCODE
1 1 Banteay Meanchey 102 Mongkol Borei 10209 Ruessei Kraok 1020914
2 1 Banteay Meanchey 104 Preah Netr Preah 10402 Chob Veari 1040201
3 1 Banteay Meanchey 107 Thma Puok 10704 Thma Puok 1070404
4 2 Battambang 202 Thma Koul 20201 Ta Pung 2020102
5 3 Kampong Cham 309 Krouch Chhmar 30905 Krouch Chhmar 3090503
6 3 Kampong Cham 310 Memot 31008 Memot 3100810
VNAME ODCODE ODNAME FACILITCOD
1 Kaoh Kaev 0101 Mongkol Borei 010001
2 Chob 0103 Preah Net Preah 010301
3 Kasen 0104 Thmar Puok 010401
4 Paoy Yong 0201 Thmar Koul 020101
5 Krouch Chhmar Leu 0304 Kroch Chhmar - Stueng Trang 030401
6 Tboung Voat 0305 Memut 030501
FACILITNAM geometry
1 Provincial Hospital POINT (286126.9 1498169)
2 Preah Net Preah POINT (302911.9 1507014)
3 Thmar Puok POINT (289389 1542390)
4 Thmar Koul POINT (293359.9 1465952)
5 Kroch Chhmar POINT (567301.9 1359213)
6 Memut POINT (628526 1308056)National Hospitals:
points_nationalhospital <- st_read(dsn = "../../data/cambodia/nationalhospital",
layer = "national_hospital_en")Reading layer `national_hospital_en' from data source
`C:\viddyasri\IS415-GAA\data\cambodia\nationalhospital' using driver `ESRI Shapefile'
Simple feature collection with 9 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 376678.3 ymin: 1276217 xmax: 491651.7 ymax: 1478968
Projected CRS: WGS 84 / UTM zone 48Nhead(points_nationalhospital)Simple feature collection with 6 features and 14 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 488043.7 ymin: 1276217 xmax: 491651.7 ymax: 1280229
Projected CRS: WGS 84 / UTM zone 48N
pcode pname dname cname vname building
1 12 Phnom Penh Daun Penh Sras Chok Not found #3
2 12 Phnom Penh Daun Penh Phsar Thmey Ti Muoy Not found #118
3 12 Phnom Penh Toul Kork Tuek L'ak Ti Pir Not found #188
4 12 Phnom Penh Chamkarmon Tumnob Tuek Not found Not found
5 12 Phnom Penh Toul Kork Tuek L'ak Ti Pir Not found #100
6 12 Phnom Penh Daun Penh Srah Chak Not found Not found
street
1 Monivong Blvd
2 Preah Norodom Blvd, corner St. Kromoun Sar
3 Yothapol Khemarak Phoumin Blvd (271)
4 Yothapol Khemarak Phoumin Blvd (271)
5 Confederation de la Russie
6 Chivapol (St.90), corner France St. (47)
web facilitcod
1 https://www.calmette.gov.kh/ NH1
2 Not found NH2
3 Not found NH3
4 http://www.khmersoviethospital.org.kh NH4
5 Not found NH5
6 https://www.beat-richner.ch/ NH6
facilitnam
1 Calmette Hospital
2 Preah Ang Duong Hospital
3 Cambodia - China Friend Preah Kossamak Hospital
4 Khmer–Soviet Friendship Hospital
5 National Pediatric Hospital
6 Kantha Bopha IV Children's Hospital
reference Lat Long language
1 screenshot_moh__17.03.2020.pdf 11.58105 104.9159 English
2 screenshot_moh__17.03.2020.pdf 11.57166 104.9234 English
3 screenshot_moh__17.03.2020.pdf 11.56397 104.8903 English
4 screenshot_moh__17.03.2020.pdf 11.54476 104.9041 English
5 screenshot_moh__17.03.2020.pdf 11.56849 104.8974 English
6 screenshot_moh__17.03.2020.pdf 11.57753 104.9213 English
geometry
1 POINT (490831.1 1280229)
2 POINT (491651.7 1279190)
3 POINT (488043.7 1278341)
4 POINT (489539.9 1276217)
5 POINT (488818 1278841)
6 POINT (491419.7 1279840)Data Preparation
Preparing Province Data:
For province data, since the province names have unique characters, we first modify them by converting characters in the NAME_1 column to ASCII equivalents using the stringi::stri_trans_general() function and putting the new, standardized names into a column called PROVINCE.
province_sf <- province_sf %>%
mutate(PROVINCE = stringi::stri_trans_general(NAME_1, "Latin-ASCII"))Preparing Health Facilities Data:
As the facilities are initially not categorised, we make a new category column and label them accordingly.
points_healthcenter <- points_healthcenter %>% mutate(CATEGORY = "Health Center")
points_healthpost <- points_healthpost %>% mutate(CATEGORY = "Health Post")
points_referralhospital <- points_referralhospital %>% mutate(CATEGORY = "Referral Hospital")
points_nationalhospital <- points_nationalhospital %>% mutate(CATEGORY = "National Hospital")Next, we remove any unwanted columns.
points_healthcenter <- subset(points_healthcenter, select = -COVERNAME)Since the columns in the national hospital dataset differ from the others, we’re starting the standardization process here. First, we’ll convert all columns to uppercase. Then, we’ll remove unmatched columns, add missing ones, and rearrange them. Finally, we’ll merge all the data points.
st_geometry(points_nationalhospital) <- "geometry"
# Get the names of all columns except the geometry column
column_names <- names(points_nationalhospital)[!grepl("^geometry$", names(points_nationalhospital))]
# Convert column names to uppercase
column_names_upper <- toupper(column_names)
# Replace column names in the sf object
names(points_nationalhospital)[!grepl("^geometry$", names(points_nationalhospital))] <- column_names_upper# Drop columns "BUILDING", "STREET", "WEB", "REFERENCE", "LAT", "LONG", and "LANGUAGE"
points_nationalhospital <- subset(points_nationalhospital, select = -c(BUILDING, STREET, WEB, REFERENCE, LAT, LONG, LANGUAGE))
points_nationalhospital$DCODE <- NA
points_nationalhospital$CCODE <- NA
points_nationalhospital$VCODE <- NA
points_nationalhospital$ODCODE <- NA
points_nationalhospital$ODNAME <- NA
# Rearrange the columns
points_nationalhospital <- points_nationalhospital[, c("PCODE", "PNAME", "DCODE", "DNAME", "CCODE", "CNAME", "VCODE", "VNAME", "ODCODE", "ODNAME", "FACILITCOD", "FACILITNAM", "CATEGORY", "geometry")]
# Print the modified sf object
print(points_nationalhospital)Simple feature collection with 9 features and 13 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 376678.3 ymin: 1276217 xmax: 491651.7 ymax: 1478968
Projected CRS: WGS 84 / UTM zone 48N
PCODE PNAME DCODE DNAME CCODE CNAME VCODE
1 12 Phnom Penh NA Daun Penh NA Sras Chok NA
2 12 Phnom Penh NA Daun Penh NA Phsar Thmey Ti Muoy NA
3 12 Phnom Penh NA Toul Kork NA Tuek L'ak Ti Pir NA
4 12 Phnom Penh NA Chamkarmon NA Tumnob Tuek NA
5 12 Phnom Penh NA Toul Kork NA Tuek L'ak Ti Pir NA
6 12 Phnom Penh NA Daun Penh NA Srah Chak NA
7 17 Siem Reap NA Siem Reap NA Kork Chak NA
8 12 Phnom Penh NA Chamkarmon NA Boeng Keng Kang 2 NA
9 12 Phnom Penh NA Daun Penh NA Wat Phnom NA
VNAME ODCODE ODNAME FACILITCOD
1 Not found NA NA NH1
2 Not found NA NA NH2
3 Not found NA NA NH3
4 Not found NA NA NH4
5 Not found NA NA NH5
6 Not found NA NA NH6
7 Trapeang Ses NA NA NH7
8 Not found NA NA NH8
9 Not found NA NA NH9
FACILITNAM
1 Calmette Hospital
2 Preah Ang Duong Hospital
3 Cambodia - China Friend Preah Kossamak Hospital
4 Khmer–Soviet Friendship Hospital
5 National Pediatric Hospital
6 Kantha Bopha IV Children's Hospital
7 Jayavarman VII Hospital
8 National Center for Tuberculosis and Leprosy Control (CENAT)
9 National Maternal and Child Health Center (NMCHC)
CATEGORY geometry
1 National Hospital POINT (490831.1 1280229)
2 National Hospital POINT (491651.7 1279190)
3 National Hospital POINT (488043.7 1278341)
4 National Hospital POINT (489539.9 1276217)
5 National Hospital POINT (488818 1278841)
6 National Hospital POINT (491419.7 1279840)
7 National Hospital POINT (376678.3 1478968)
8 National Hospital POINT (491266.3 1277206)
9 National Hospital POINT (491327.1 1280024)points_facilities <- rbind(points_healthcenter, points_healthpost, points_referralhospital, points_nationalhospital)points_facilities$PNAME <- gsub("Banteay Mean Chey", "Banteay Meanchey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Battambang", "Batdambang", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Battam Bang", "Batdambang", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Speu", "Kampong Spoe", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Spueu", "Kampong Spoe", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kampong Thom", "Kampong Thum", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Koh Kong", "Kaoh Kong", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Kratie", "Kracheh", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Pailin", "Krong Pailin", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Sihanoukville", "Krong Preah Sihanouk", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Sihaknouk Vill", "Krong Preah Sihanouk", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Mondul Kiri", "Mondol Kiri", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddor Meanchey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddar Mean chey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Oddar Meanchey", "Otdar Mean Chey", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Pursat", "Pouthisat", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Ratanak Kiri", "Rotanokiri", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Siemreap", "Siemreab", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Siem Reap", "Siemreab", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Stung Treng", "Stoeng Treng", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Stung Treng", "Stoeng Treng", points_facilities$PNAME)
points_facilities$PNAME <- gsub("Takeo", "Takev", points_facilities$PNAME)st_geometry(points_facilities)Geometry set for 1130 features
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 213103.9 ymin: 1155396 xmax: 764437 ymax: 1594142
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:write_rds(points_facilities, "../../data/rds/points_facilities.rds")Handle Invalid Geometries:
length(which(st_is_valid(province_sf) == FALSE))[1] 0length(which(st_is_valid(cam_road_sf) == FALSE))[1] 0length(which(st_is_valid(points_facilities) == FALSE))[1] 0Projection System Transformation:
province_sf <- st_transform(province_sf, 32648)
st_crs(province_sf)Coordinate Reference System:
User input: EPSG:32648
wkt:
PROJCRS["WGS 84 / UTM zone 48N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 48N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",105,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 102°E and 108°E, northern hemisphere between equator and 84°N, onshore and offshore. Cambodia. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Russian Federation. Singapore. Thailand. Vietnam."],
BBOX[0,102,84,108]],
ID["EPSG",32648]]st_geometry(province_sf)Geometry set for 25 features
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 211598 ymin: 1096586 xmax: 784872.1 ymax: 1625374
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:write_rds(province_sf, "../../data/rds/province_sf.rds")cam_road_sf <- st_transform(cam_road_sf, 32648)
st_crs(cam_road_sf)Coordinate Reference System:
User input: EPSG:32648
wkt:
PROJCRS["WGS 84 / UTM zone 48N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 48N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",105,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 102°E and 108°E, northern hemisphere between equator and 84°N, onshore and offshore. Cambodia. China. Indonesia. Laos. Malaysia - West Malaysia. Mongolia. Russian Federation. Singapore. Thailand. Vietnam."],
BBOX[0,102,84,108]],
ID["EPSG",32648]]st_geometry(cam_road_sf)Geometry set for 169682 features
Geometry type: MULTILINESTRING
Dimension: XY
Bounding box: xmin: 212247.3 ymin: 1152259 xmax: 783030 ymax: 1625600
Projected CRS: WGS 84 / UTM zone 48N
First 5 geometries:write_rds(cam_road_sf, "../../data/rds/cam_road_sf.rds")Network Constrained Kernel Density Estimation
For the Shiny application, all kernel types will be included and the user will be able to experiment with different fixed bandwidths. NetKDE can also be performed on the 25 provinces of Cambodia and can be narrowed down to facility type. However, for the sake of execution time, this exercise will demonstrate only 3 combinations of kernels and bandwidths for the simple method specifically for the province of Kampot.
Preparing the Road Layer and Facility Point Data:
province_sf <- read_rds("../../data/rds/province_sf.rds")
points_facilities <- read_rds("../../data/rds/points_facilities.rds")
cam_road_sf <- read_rds("../../data/rds/cam_road_sf.rds")province = province_sf %>% filter(PROVINCE=="Kampot")
province_road_sf = st_intersection(cam_road_sf, st_union(province))province_facilities = st_intersection(points_facilities, st_union(province))Converting to Simple Geometries:
province_road_sf <- st_cast(province_road_sf, "LINESTRING")Basic Plot:
tmap_mode("plot")
tm_shape(province_road_sf) +
tm_lines(col = "#2b2b2b") +
tm_shape(province_facilities) +
tm_dots(col = "#29ab87", size = 0.2)
Performing NetKDE
Prepare Lixel Objects and Line Centre Points:
lixels <- lixelize_lines(province_road_sf,
750,
mindist = 375)
samples <- lines_center(lixels)Quartic Kernel
Bandwidth - 200:
densities_q200 <- nkde(
province_road_sf,
events = province_facilities,
w = rep(1, nrow(province_facilities)),
samples = samples,
kernel_name = "quartic",
bw = 200,
div = "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1, 1),
max_depth = 8,
agg = 5,
sparse = TRUE,
verbose = FALSE
)samples$density_q200 <- densities_q200*nrow(province_facilities)*1000
lixels$density_q200 <- densities_q200*nrow(province_facilities)*1000Visualization:
samples2 <- samples[order(samples$density_q200),]
tmap_mode("plot")
tm_shape(province_road_sf) +
tm_lines("black") +
tm_shape(samples2) +
tm_dots("density_q200", style = "kmeans", palette = "GnBu", n = 7, size = 0.07) +
tm_layout(legend.outside = FALSE)
province_facilities <- province_facilities %>%
mutate(COMBINED_ID = paste("Name:", FACILITNAM, "|| Category:", CATEGORY))
tmap_mode('view')
tm_basemap(server = "Esri.WorldTopoMap") +
tm_basemap(server = "Esri.WorldGrayCanvas")+
tm_basemap(server = "OpenStreetMap") +
tm_shape(lixels) +
tm_lines(col="density_q200", palette="PuRd", lwd=5) +
tm_shape(province_facilities) +
tm_dots(id = "COMBINED_ID")+
tm_shape(province) +
tm_borders()